\bfi % ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff

\mbox{%m
\hspace{-5ex}
a)\raisebox{-0.35\textwidth}{\ig{height=0.35\textwidth}{lpconv.eps}}
b)\raisebox{-0.35\textwidth}{\ig{height=0.35\textwidth}{lpeig.eps}}
c)\raisebox{-0.35\textwidth}{\ig{height=0.35\textwidth}{lpeig_bwlp.eps}}
}%m
\ca{a) Pointwise convergence of exterior Helmholtz Dirichlet BVP with
double-layer representation, at $k=10$ (solid: pure DLP, dotted,
combined-field formulation).
b) Spectrum of $(\frac{1}{2}I + D)$ in the complex plane (for $M=80$).
We have superimposed a circle of radius 1/2 centered at 1/2, to which the
eigenvalues appear to
condense in the $k\to\infty$ limit---to prove this is an open problem!
c) Spectrum of $(\frac{1}{2}I + D - ikS)$ relevant for combined-field
formulation.
}{f:lp}
\efi

% ----------------------------------------------------------------------------
\section{Layer potentials}
\label{s:layer}

Layer potentials are representations of Helmholtz solutions
involving a weighted integral of a fundamental solution
over a boundary \cite{coltonkress}.
They are similar to the MFS representations discussed above, with the
crucial advantage that a {\em second-kind}
formulation is often possible, i.e. the operator involved is
identity plus compact, and the resulting matrix problems are well-conditioned.
It is easy to set up single- and double-layer potentials (SLP and
DLP) living on a segment.
Starting with the trefoil segment {\tt tref} defined in Sec.~\ref{s:ext}, a DLP
representation in its exterior is set up by,
\begin{verbatim}
tref = segment.radialfunc([], {@(q) 1 + 0.3*cos(3*q), @(q) -0.9*sin(3*q),...
                   @(q) -2.7*cos(3*q)});   % note R" also given, helps later
d = domain([], [], tref, -1); d.k = 10;    % exterior domain & wavenumber
d.addlayerpot(tref, 'D');                  % adds DLP on tref segment
\end{verbatim}
The coefficients in {\tt p.co} represent {\em values} of the double-layer
density function at the quadrature points.
In contrast to MFS bases,
the segment's own quadrature points are always used. Therefore the
number of degrees of freedom is the current value of {\tt numel(tref.x)},
and when you update $N$ for a {\tt layerpot} you actually
change the quadrature point number of the underlying segment.

We may repeat the exterior BVP convergence of Sec.~\ref{s:ext} using the DLP,
\begin{verbatim}
f = @(z) besselh(0, d.k * abs(z-0.3-0.2i));  % known exterior field
tref.setbc(1, 'D', [], @(t) f(tref.Z(t)));   % its Dirichlet data
p = bvp(d);
Ns = 5:5:80; for i=1:numel(Ns)
  p.updateN(Ns(i)); p.solvecoeffs; N(i) = p.N;
  e(i) = f(2) - p.pointsolution(pointset(2));   % error at a point in d
end
figure; semilogy(N, abs(e), '+-'); xlabel('N'); ylabel('error in u(2)');
\end{verbatim}
The result is the solid line in Fig.~\ref{f:lp}a.
The condition number of the problem is $O(1)$, viz.\ {\tt cond(p.A)} gives
69, as opposed to the ill-conditioned value
$1.8\times 10^7$ for MFS in Sec.~\ref{s:ext}.
Note that errors can no longer be estimated by {\tt p.bcresidualnorm} since it
always around machine precision, $10^{-16}$.

What actually happened above? Quite a lot! The BVP class needed to know how
to evaluate the DLP on its own source segment, which requires two new
ingredients:
\ben
\item the {\em jump relation} that states that
$u^+ = (\frac{1}{2}I + D)\tau$ where $\tau$ is the density function and
$D$ the integral operator on the boundary with the double-layer kernel, and,
\item a good quadrature scheme to evaluate $D\tau$ at the quadrature points
themselves is needed; by default for periodic quadrature%
  \footnote{\mpspack\ also implements Kapur-Rokhlin schemes of order 2, 4,
    and 10, which for non-FMM-accelerated problems sich as this are inferior
    to Kress'. For non-periodic open segments, no special quadratures are
    currently implemented.}
the Martensen--Kussmaul scheme recommended by Kress is used \cite{coltonkress},
and this requires $z''(t)$, hence the need for $R''$ in {\tt radialfunc}.
\een
The above scheme is called the {\em Nystr\"{o}m method} \cite{coltonkress}.
How did \mpspack\ know to evaluate $u^+$ rather than $u^-$? Only the
$+$ side of the segment is connected to the evaluation domain.
Note that our choice of normalization is that the matrix
{\tt p.A} maps density values to field values multiplied by the square-root
of quadrature weights. To recover the unweighted matrix and check its
eigenvalues, which approximate well eigenvalues of $(\frac{1}{2}I + D)$,
we use
\co{figure; plot(eig(diag(1./p.sqrtwei)*p.A), '+'); axis equal}
giving Fig.~\ref{f:lp}b.

This plot shows an eigenvalue alarmingly close to zero. Why?
It is well-known that $\frac{1}{2}I + D$ goes singular when $k^2$ is
a Dirichlet eigenvalue of the interior domain \cite{coltonkress}.
This leads to blow-up of the condition number, hence non-robustness.
A cure due to Brakhage, Werner, Leis and Panich \cite{coltonkress}
is to use a combined-field representation of DLP minus $ik$ times SLP,
which is done by replacing the above basis setup by
\co{d.addlayerpot(tref, [-1i*d.k 1]);}
Delightfully, the spectrum of $\frac{1}{2}I + D - ikS$ avoids
zero like the plague (Fig.~\ref{f:lp}c) and now {\tt cond(p.A)} is reduced to 5.

\subsection{Transmission problems}
\label{s:trans}

We return now to the transmission problem from the end of Section~\ref{s:scatt}.
The Rokhlin--M\"{u}ller scheme \cite{rokh83}
for scattering from a dielectric obstacle
uses independent single- and double-layer densities on the obstacle
boundary curve, each of which is used to represent the field in
both interior and exterior.
This means that a single density must affect
{\em two} domains in which, however, the wavenumbers are different.
This is easy to set-up in \mpspack, and it contrasts the cases considered
until now where a basis object may affect only a single domain.
The advantage of this scheme is that
the hypersingular kernels on the plus and minus sides of the surface
cancel to give a second-kind Fredholm system.
Some code for the below is found in {\tt examples/dielscatrokh.m}

We create interior and exterior domains from a single segment
(in fact the same one as in Sec.~\ref{s:ext}, using a simpler
command), and
enforce continuity conditions at their interface, as follows,
\begin{verbatim}
M = 330; s = segment.smoothstar(M, 0.3, 3);       % smooth closed segment
di = domain(s, 1); di.setrefractiveindex(n);      % interior
de = domain([], [], s, -1);                       % exterior
setmatch(s, 'diel', 'TM');
\end{verbatim}
We use a new command to add `double-sided' layer-potentials to the segment,
one DLP and one SLP,
\begin{verbatim}
s.addinoutlayerpots('d'); s.addinoutlayerpots('s');
\end{verbatim}
We may check that the DLP \verb+di.bas{1}+ affects both domains by
checking that \verb+di.bas{1}.doms+ is a list of two domains.
By default, a Kress quadrature scheme is used for closed segments---this
can be changed in the options to the above calls.%
  \footnote{currently, since they cancel,
    the hypersingular parts of the $T$ operator
    are not even included in the $T$ evaluation}.
Given the above domains (and once boundary/matching conditions and
basis sets have been chosen), we set up then solve the problem
with the same parameters
as in Section~\ref{s:scatt},
\begin{verbatim}
p = scattering(de, di);
p.setoverallwavenumber(30);
p.setincidentwave(pi/6);              % if just angle given, it's a plane wave
p.solvecoeffs;
\end{verbatim}
A call to {\tt p.showthreefields} now reproduces Fig.~\ref{f:diel}.
Note that in order to get 13 digits of accuracy at this medium wavenumber,
we needed $M=330$, i.e.\ 660 degrees of freedom. (This is larger than the number
needed in Section~\ref{s:scatt}, although the accuracy was lower there).
At $k=8$, $M=110$ is sufficient for 14 digit accuracy.
However, the linear system condition number {\tt cond(p.A)} is only
$6\times 10^3$, enabling the use of iterative methods such as GMRES,
in contrast to Sec.~\ref{s:scatt} in which it
was highly ill-conditioned. The convergence rate is exponential.


In a future edition we will discuss
sound-hard scattering and direct formulations.
Please get in touch with the authors if you have any questions
or suggestions!
